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1. Introduction 

Understanding the dynamical origin of the mechanisms which underly the phenomenol- 
ogy of heat conduction has remained one of the major open problems of statistical 
mechanics ever since Fourier's seminal work pQ. Fourier himself actually warned his 
reader that the effects of heat conduction "make up a special range of phenomena 
which cannot be explained by the principles of motion and equilibrium," thus seemingly 
rejecting the possibility of a fundamental level of description. Of course, with the sub- 
sequent developments of the molecular kinetic theory of heat, starting with the works 
of the founding-fathers of statistical mechanics, Boltzmann, Gibbs and Maxwell, and 
later with the definite triumph of atomism, thanks to Perrin's 1908 measurement of the 
Avogadro number in relation to Einstein's work on Brownian motion, Fourier's earlier 
perception was soon discarded as it became clear that heat transport was indeed the 
effect of mechanical causes. 

Yet, after well over a century of hard labor, the community is still actively hunting 
for a first principles based derivation of Fourier's law, some authors going so far as to 
promise "a bottle of very good wine to anyone who provides [a satisfactory answer to 
this challenge]" [2J. Thus the challenge is, starting from the Hamiltonian time evolution 
of a system of interacting particles which models a fluid or a crystal, to derive the 
conditions under which the heat flux and temperature gradients are linearly related by 
the coefficient of heat conductivity. This is the embodiment of Fourier's law. 

As outlined in [2] , it is necessary, in order to achieve this, that the model statistical 
properties be fully determined in terms of the local temperature, a notion which involves 
that of local thermalization. To visualize this, we imagine that the system is divided 
up into a large number of small volume elements, each large enough to contain a 
number of particles whose statistics is accurately described by the equilibrium statistics 
at the local temperature associated to the volume element under consideration. To 
establish this property, one must ensure that time scales separate, which implies that 
the volume elements settle to local thermal equilibrium on time scales much shorter 
than the ones which characterize the transport of heat at the macroscopic scale. Local 
thermal equilibrium therefore relies on very strong ergodic properties of the model. 

The natural framework to apply this programme is that of chaotic billiards. Indeed 
non-interacting particle billiards, where individual tracers move and make specular 
collisions among a periodic array of fixed convex scatterers (the periodic Lorentz gases) 
are the only known Hamiltonian systems for which mass transport [3] and shear and bulk 
viscosities [4] have been established in a rigorous way. Here, rather than local thermal 
equilibrium, it is a relaxation to local equilibrium, occurring on the constant energy sheet 
of the individual tracer particles, which allows to model the mass transport by a random 
walk and yields an analytic estimate of the diffusion coefficient [5]. Arguably Lorentz 
gases, whether periodic or disordered, in or out of equilibrium, have played, over more 
than a century, a privileged and most important role in the development of transport and 
kinetic theories [6]. However, in the absence of interaction among the tracer particles, 
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there is no mechanism for energy exchange and therefore no process of thermalization 
before heat is conducted. If, on the other hand, one adds some interaction between 
the tracers, say, as suggested in [2], by assigning them a size, the resulting system will 
typically have transport equations for the diffusion of both mass and energy. Though 
these systems may be conceptually simple, they are usually mathematically too difficult 
to handle with the appropriate level or rigor. 

An example of billiard with interacting tracer particles is the modified Lorentz gas 
with rotating discs proposed by Mejia-Monasterio et al. [7j. This system exhibits normal 
transport with non-trivial coupled equations for heat and mass transport. Though this 
model has attracted much attention among mathematicians [El El [TUl EEI1 E21 EE], the 
rigorous proof of Fourier's law for such a system of arbitrary size and number of tracers 
relies on assumptions which, though they are plausible, are themselves not resolved. 
Moreover the equations which govern the collisions and the discs rotations do not, as 
far as we know, derive from a Hamiltonian description. 

Yet a remedy to such limitations with respect to the number of particles in 
interaction that a rigorous treatment will handle may have been just around the corner 
[14] , and, this, within a Hamiltonian framework. Indeed in [15], Bunimovich et al. 
introduced a class of dispersing billiard tables with particles that are in geometric 
confinement —i. e. trapped within cells- but that can nevertheless interact among 
particles belonging to neighbouring cells. The authors proved ergodicity and strong 
chaotic properties of such systems with arbitrary number of particles. 

More recently, the idea that energy transport can be modeled as a slow diffusion 
process resulting from the coupling of fast energy- conserving dynamics has led to proofs 
of central limit theorems in the context of models of random walks and coupled maps 
which describe the diffusion of energy in a strongly chaotic, fast changing environment 
[T6| [T7] . Although the extension of these results to symplectic coupled maps, let alone 
Hamiltonian flows, is not yet on the horizon, it is our belief that such systems as the 
dispersing billiard tables introduced in [15] will, if any, lend themselves to a fully rigorous 
treatment of heat transport within a Hamiltonian framework. As we announced in 
[18] . the reason why we should be so hopeful is that the particles confinement has two 
important consequences : first, relaxation to local thermal equilibrium is preceded by a 
relaxation of individual particles to local equilibrium^, which occurs at constant energy 
within each cell, and has strong ergodic properties that guarantee the rapid decay of 
statistical correlations; and, second, heat transport, unlike in the rotating discs model, 
can be controlled by the mere geometry of the billiard, which also controls the absence 
of mass transport. As of the first property, the relaxation to a local equilibrium before 
energy exchanges take place is characterized by a fast time scale, much faster than that 
of relaxation to local thermal equilibrium among neighbouring cells, which is itself much 
faster than the hydrodynamic relaxation scale. There is therefore a hierarchy of three 

% Let us underline the distinction we make here and in the sequel between relaxation to local 
equilibrium, which precedes energy exchanges, and relaxation to local thermal equilibrium, which 
involves energy exchanges among particles belonging to neighbouring cells. 
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separate time scales in this system, the first accounting for relaxation at the microscopic 
scale of individual cells, the second one at the mesoscopic scale of neighbouring cells, 
and the third one at the macroscopic scale of the whole system. 

In this paper, we achieve two important milestones towards a complete first 
principles derivation of the transport properties of such models. Having defined the 
model, we establish the conditions for separation of time scales and relaxation to local 
equilibrium, identifying a critical geometry where binary collisions become impossible. 
Assuming relaxation to local equilibrium holds, we go on to considering the time 
evolution of phase-space densities and derive, from it, a master equation which governs 
the exchange of energy in the system [19], thus going from a microscopic scale description 
of the Liouville equation to the mesoscopic scale at which energy transport takes place. 
We regard this as the first milestone, namely identifying the conditions under which 
one can rigorously reduce the level of descrition from the deterministic dynamics at 
the microscopic level to a stochastic process described by a master equation at the 
mesoscopic level of energy exchanges. 

This master equation is then used to compute the frequency of binary collisions 
and to derive Fourier's law and the macroscopic heat equation, which results from the 
application of a small temperature gradient between neighbouring cells. This is our 
second milestone : an analytic formula for heat conductivity, exact for the stochastic 
system, and thus valid for the determinisitc system at the critical geometry limit. 
These results are then checked against numerical computations of these quantities, with 
outstanding agreement, and shown to extend beyond the critical geometry, with very 
good accuracy, to a wide range of parameter values. 

We further characterize the chaotic properties of the model and offer arguments 
to account for the spectrum of Lyapunov exponents of the system, as well as the 
Kolmogorov-Sinai entropy, expressions which are exact at the critical geometry. Again, 
these results are very nicely confirmed by our numerical computations. 

The paper is organized as follows. The models, which we coin lattice billiards, are 
introduced in section [2j Their main geometric properties are established, distinguishing 
transitions between insulating and conducting regimes under the tuning of a single 
parameter. The same parameter controls the time scales separation responsible for local 
equilibrium. Section [3] provides the derivation of the master equation which, under the 
assumption of local equilibrium, governs energy transport. The main observables are 
computed and their scaling properties discussed. In section H] we review the properties 
of the model and assess the validity of the results of section [3] under the scope of our 
numerical computations. The chaotic properties of the models are discussed in section 
EJ We use simple theoretical arguments to predict some of these properties and compare 
them to numerical computations. Finally, conclusions are drawn in section [6j 
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2. Lattice billiards 

To introduce our model, we start by considering the uniform motion of a point particle 
about a dispersing billiard table, B p , defined by the domain exterior to four overlapping 
discs of radii p, centered at the four corners of a square of sides I. The radius is thus 
restricted to the interval 1/2 < p < l/y/2, where the lower bound is the overlap (or 
bounded horizon) condition, and the upper bound is reached when B p is empty. 



Figure 1. Two equivalent representations of a dispersing billiard table : (left) a point 
particle moves uniformly inside the domain B p and performs specular collisions with 
its boundary; (right) a disc of radius p m moves uniformly inside the domain B p and 
performs specular collisions with fixed discs of radii pf = p — p m . The radius p m is a 
free parameter which is allowed to take any value between and p. 



As illustrated in figure dj the motion of a point particle in this environment is 
a limiting case of a class of equivalent dispersing billiards, whereby the point particle 
becomes a moving disc with radius p m , < p m < p, and bounces off fixed discs of radii 
Pf = p — p m . In all these cases, the motion of the center of the moving disc is equivalent 
to that of the point particle in B p . The border of the domain B p constitutes the walls 
of the billiards. 

In the absence of cusps (which occur at p = 1/2), the ergodic and hyperbolic 
properties of these billiards are well established [20]. In particular, the long term 
statistics of the billiard map, which takes the particle from one collision event to the 
next, preserves the measure cos0drd0, where <fi denotes the angle that the particle 
post-collisional velocity makes with respect to the normal vector to the boundary. A 
direct consequence of this invariance is a general formula which relates the mean free 
path, £, to the billiard table area, \B P \, and perimeter \dB p \, £ = ir\B p \/\dB p \. 

The ratio between the speed of the particle, which we denote v, and the mean free 
path gives the wall collision frequency[§], v c = v/ £, whose computation is shown in figure 
[2j With this quantity, one can relate the billiard map iterations to the time- continuous 
dynamics of the flow. In particular, the billiard map has two Lyapunov exponents, 
opposite in signs and equal in magnitudes, which, multiplied by the wall collision 

§ Anticipating the more general definition of the wall collision frequency for interacting particle billiard 
cells, we adopt the subscript "c" in reference to "critical" for reasons to be clarified below. 
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frequency, correspond to the two non-zero Lyapunov exponents of the flow (there are 
two additional zero exponents related to the direction of the flow and conservation of 
energy). The results of numerical computations of the positive one, denoted A+, are 
shown in figure [2] for different values of the parameters p. 
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Figure 2. Collision frequency v c (solid line) and numerical computation of the 
corresponding positive Lyapunov exponent of the flow (dots) of dispersing billiards 
such as shown in figure [TJ Here we took v = 1 and the square side to be I = l/y/2. 

Thus let Q(i,j) denote the rhombus of sides I centered at point 

(y/2li,lj/y/2), j even, 

(V2U + l/y/2, lj/V2), 3 odd. 1 ' 

The rhombic billiard cell illustrated in figure [T] becomes a domain centered at point 
(cij,dij), defined according to 

Bpihj) = \( x ,v) e Q{hj) 8[(x,y),(c ij +p k ,d ij +q k )} > p,k = l,...,d},(2) 



[piji dij ) 



where S[., .) denotes the usual Euclidean distance between two points, and {pk,Qk) = 
(±//v^2, 0), (0, ±Z/a/2) are the coordinates of the d = 4 discs at the corners of the 
rhombus. 

Now consider a number of copies {B p (i, which tessellate a two-dimensional 
domain. We define a lattice billiard as a collection of billiard cells 



(3) 



VI <i,i' <m,l <j,j'<n 2 ,i^i',j^j'}. 



Each individual cell of this billiard table possesses a single moving particle of radius p m , 
< p m < P and unit mass. All the moving particles are assumed to have independent 
initial coordinates within their respective cells, with the proviso that no overlap can 
occur between any pair of moving particles. The system energy, E = Nk^T (with 
N = ni x n 2 , the number of moving particles, T the system temperature, and k-Q 
Boltzmann's constant) is constant and assumed to be initially randomly divided among 
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the kinetic energies of the moving particles, E = £\ . ey, ey = mvfj/2, where %■ denotes 
the speed of particle 

Energy exchanges occur when two moving particles located in neighbouring cells 
collide. Such events can take place provided the radii of the moving particles p m is 
large enough compared to p. Indeed the value of the critical radius, below which binary 
collisions do not occur, is determined by half the separation between the corners of two 
neighbouring cells, 




Figure 3. A binary collision event in the critical configuration where p m = p c would 
occur only provided the colliding particles visit the corresponding corners of their cells 
simultaneously. The value of p in this figure is the same as in figures [T] and 2J 

For the sake of illustration, the unlikely occurrence of a binary collision event at 
the critical radius p m = p c is shown in figure [3j 

All the collisions are elastic and conserve energy, so that the dynamical system is 
Hamiltonian with 2N degrees of freedom. Its phase space of positions and velocities 
is 4iV-dimensional. Accordingly, the sensitivity to initial conditions of the dynamics is 
characterized by 4iV Lyapunov exponents, {A,}^, obeying the pairing rule of symplectic 
systems, A4Ar_ i+ i = —A,,, i — 1, . . . , 2N. 

Collision events between two moving particles are referred to as binary collision 
events and will be distinguished from wall collision events, which occur between the 
moving particles and the walls of their respective confining cells. The occurrences of 
the former are characterized by a binary collision frequency, v^, and the latter by a wall 
collision frequency, v w . Both frequencies depend on the difference p m — p c , separating 
the moving particles radii from the critical radius, equation (Tj0). By definition of p c , 
the binary collision frequency vanishes at p m = p c , u h \ Pm=Pc = 0, and, correspondingly, 
the wall collision frequency at the critical radius is the collision frequency of the single- 
cell billiard, v w \p m =p c — v c . We will assume from now, unless otherwise stated, that 
the system is globally isolated and apply periodic boundary conditions at the borders, 
thereby identifying B p {i + kni,j + ln 2 ) with B p (i,j) for any k, I £ Z, 1 < i < n\, 
1 < j < n 2 . 

Examples of such billiards are displayed in figure HI Obviously the quincunx 
rhombic lattice structure, which is generated by the rhombic cells, is but one among 
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Figure 4. Examples of lattice billiards with triangular (top), rhombic (middle) and 
hexagonal (bottom) tilings. The coloured particles move among an array of fixed 
black discs. The radii of both fixed and mobile discs are chosen so that (i) every 
moving particle is geometrically confined to its own billiard cell (identified as the area 
delimited by the exterior intersection of the black circles around the fixed discs), but 
(ii) can nevertheless exchange energy with the moving particles in the neighbouring 
cells through binary collisions. The solid broken lines show the trajectories of the 
moving particles centers about their respective cells. The colours are coded according 
to the particles kinetic temperatures (from blue to red with increasing temperature). 
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different possible structures. Triangular, upright square, or hexagonal cells can be used 
as alternative periodic structures. One might also cover the plane with random or 
quasi- crystalline tessellations. The only relevant assumptions in what follows is that 
the moving particles must be confined to their (dispersing) billiard cells and that binary 
interactions between neighbouring cells can be turned on and off by tuning the system 
parameters. 

The two important features of such lattice billiards is that (i) there is no mass 
transport across the billiard cells since the moving particles are confined to their 
respective cells, and (ii) energy transport can occur through binary collision events 
which take place when the particles of two neighbouring cells come into contact. In 
periodic structures such as the quincunx rhombic lattice, the possibility of such collisions 
is controlled by tuning the parameter p m about the critical radius p c , keeping p fixed. 

We can therefore distinguish two separate regimes : 

• Insulating billiard cells: < p m < p c 

Absence of interaction between the moving discs. No transport process across the 
individual cells can happen; 

• Conducting billiard cells: p c < p m < p 

Binary collision events are possible. Energy transport across the individual cells 
takes place. 

The case p m = p c is singular. We will refer to the critical geometry as the limit p m — > p c . 

In the insulating regime, there is no interaction among moving particles so that 
the billiard cells are decoupled. The moving particles are independent and their kinetic 
energies are individually conserved, resulting in 2N zero Lyapunov exponents. The 
equilibrium measure in turn has a product structure and phase-space distributions are 
locally uniform with respect to the particles positions and velocity directions. The 
N positive Lyapunov exponents of the system are all equal to the positive Lyapunov 
exponent of the single-cell dispersing billiard, up to a factor corresponding to the 
particles speeds, Vij = \j2eij/m : Ay, = Vij\ + , where A + is the Lyapunov exponent of 
the single-cell billiard measured per unit length. 

When particles are allowed to interact, on the other hand, local energies are 
exchanged through collision events. Thus only the total energy is conserved in the 
conducting regime. The ergodicity of such systems of geometrically confined particles 
in interaction was proven by Bunimovich et al. [15J. The resulting dynamical system, 
whose equilibrium measure is the microcanonical one (taking into consideration that 
particles are otherwise uniformly distributed within their respective cells), enjoys the 
X-property. This implies ergodicity, mixing, and strong chaotic properties, including 
the positivity of the Kolmogorov- Sinai entropy. Two Lyapunov exponents are zero, one 
associated to the conservation of energy, the other to the direction of the flow. The 
2(2N — 1) remaining Lyapunov exponents form non- vanishing pairs of exponents with 
opposite signs, Ai > . . . > \2N-1 > 0, A4Ar_i + i = — Aj, i = 1, . . . , 2N — 1. 

The regime of interest to us is that corresponding to particles interacting rarely, 
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which is to say, in analogy with a solid, that particles mostly vibrate inside their cells, 
ignorant of each other, and only seldom making collisions with their neighbours, thereby 
exchanging energy. As we turn on the interaction and let p m > p c , binary collisions, 
though they can occur, will remain unlikely. This is to say that the binary collision 
frequency, u^, will, in this regime, remain small with respect to the wall collision 
frequency, u w , which in the absence of interaction and, in particular, at the critical 
geometry, we recall is equal to the wall collision frequency of the single-cell billiard, 
v w \ Pm = Pc = v c . When p m > p c , we therefore expect z/ w ^> z/b, as well as u w ~ u c . In 
words : time scales separate. The consequence is that relaxation to local equilibrium 
-i. e. uniformization of the distribution of the particles positions and velocity directions 
at fixed speeds- occurs typically much faster than the energy exchange which drives the 
relaxation to the global equilibrium. This mechanism justifies resorting to kinetic theory 
in order to compute the transport properties of the model. 

3. Kinetic theory 

3.1. From Liouville's equation to the master equation 

The phase-space probability density is specified by the iV-particles distribution function 
p N (ri,Vx, . . . ,r N ,v N ,t), where r a and v a , a = 1, . . . , N, denote the ath particle position 
and velocity vectors. The index a stands for the label of the cells defined by 

equation (T5]). For our system, as is customary for hard sphere dynamics, this distribution 
satisfies a pseudo-Liouville equation [21], which is well defined despite the singularity 
of the hard-core interactions. This equation, which describes the time evolution of p^ 
is composed of three types of terms: (i) the advection terms, which account for the 
displacement of the moving particles within their respective billiard cells; (ii) the wall 
collision terms, which account for the wall collision events, between the moving particles 
and the d fixed scattering discs which form the cells walls; and (iii) the binary collision 
terms, which account for binary collision events, between moving particles belonging to 
neighbouring billiard cells : 



Each wall collision term involves a single moving particle with index a and one of the d 
fixed discs in the corresponding cell, with index k and position R k . Let r ak = r a — R k 
denote their relative position. Following [22], we have 




where e denotes the normal unit vector to the fixed disc k in the cell of particle a. 







(6) 
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Likewise the binary collision operator, written in terms of the relative positions r ab 
and velocities v ab of particles a and b, and the unit vector e ab that connects them, is 

B M p N (. ..,r a ,v a ,..., n, v b , ■■■) = 



/ de ab (e ab ■ v ab ) S(r ab - 2p m e ab 

J e„h V„ h >0 



2p r 

x Pn{ ...,r a ,v a - e ab (e ab ■ v ab ), ...,r b ,v b + e ab (e ab ■ v ab ), . . . ) 

- 5(r ab + 2p ra e ab )p N ( . . . , r a , v a , . . . , r b , v b , . . . ) . (7) 

We notice that only the terms B^ a ' b ^ corresponding to first neighbours are non- vanishing 
and contribute to the double sum on the RHS of equation (jSJ). 

Provided we have a separation of time scales between wall and binary collisions, the 
advection and wall collision terms on the RHS of equation ((Sj) will typically dominate 
the dynamics on the short time ~ l/v w , which follows every binary collision event, thus 
ensuring, thanks to the mixing within individual billiard cells, the relaxation of the 
phase-space distribution p^ to local equilibrium well before the occurrence of the next 
binary event, whose time scale is ~ 1/Vb- In other words, Piv(^i, Vi, ■ ■ ■ , fN, V N, t) 
quickly relaxes to a locally uniform distribution, which depends only on the local 
energies, justifying the introduction of 

f N N 

P^ eq) (ei, . . . ,e N ,t) = / Y[d r adv a p N (r 1 ,v 1 , . . . ,r N ,v N ,t)Y[5(e a - mvl/2), (8) 

o=l a=l 

where v a = |v |. On the time scale of binary collision events, this distribution 
subsequently relaxes to the global microcanonical equilibrium distribution. This process 
accounts for the transport of energy, and can be characterized by the master equation 

m 

N f 



x 



a,6=l ' 

,(leq), 



W(e a + n,e b - r]\e a ,e b )Px >(■ ■ ■ , e a + rj, . . . , e b - r], . . . , t) 

-W(e a , e b \e a -r],e b + r/)P^ eq) (. . . , e OJ . . . , e b , . . . , t) , (9) 

where W(e a ,e b \e a — rj, e b + r]) denotes the probability that an energy r] be transferred 

from particle a to particle b as the result of a binary collision event between them. 

This equation is a closure for the local equilibrium distribution P^ eq \ obtained from 

equation (jSJ) under the assumption that u w ^> u b . The first two terms on the RHS 

of equation (jSJ) are eliminated because they leave invariant the local distribution 

0^1 $( e a ~~ mv a/2)- There remain the contributions (J7|) from the binary collisions, 

which, under the assumption that the local distibutions are uniform with respect to the 

positions and velocity directions, yield the following expression of W : 

2p m m 2 f f 

W(e a ,e b \e a -r],e b + T]) = " — ry / d<pdR / dv a dv b (10) 

(2vr) 2 1 £ PiPm (2) | J Je ab . Vab >o 

x e ab ■ v ab 6 (e a - ^t»^J 5 (e b - ^vfj 6 (j) - y \(e ab ■ v a ) 2 - (e ab ■ v b ) 2 ]^ , 
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where the first integration is performed over the positions of the center of mass, 
R = (r a + r&)/2, between the two particles a and b, given that they are in contact and 
both located in their respective cells, and over the angle of the unit vector connecting 
a and b, e ab = (cos 0, sin 0). The normalizing factor \£ PtPm (2)\ denotes the 4- volume of 
the billiard corresponding to two neighbouring cells a and b, which, with the assumption 
that p m > p c , can be approximated by |£ PjPm (2)| ~ \B P \ 2 . This substitution amounts 
to neglecting the overlap between the two particles; see equation (123]) for a refinement 
of that approximation. We point out that the position and velocity integrations in 
equation (jTUI) can be formally decoupled; in this way, we can prove that the transition 



rate W is given in terms of Jacobian elliptic functions, see Appendix A 



3.2. Geometric factor 

As we show in Appendix A , an important property of the master equation (Q is that 



the factor which accounts for the geometry of collision events factorizes from the part of 
the kernel that accounts for energy exchanges. Therefore, as p m —>■ p c , the critical value 
of the radius at which binary collision events become impossible, which is the regime 
where the billiard properties are accurately described in terms of the master equation 
above, the geometric factor J d<pdR encloses the scaling properties of observables with 
respect to the billiard geometry. We now compute this quantity. 

A binary collision occurs when particles a and b come to a distance 2p m of each 
other, with r a G B p (a) and r& G B p {b). Let R = (x,y) be the center of mass coordinates 
and be the angle between the particles relative position and the axis connecting the 
center of the cells. Taking a reference frame centered between the cells, we may write 

ri = -(x,y) + o-ip m (cos(f), sin 0), (11) 

where Oi = ±1 and i = a or b. The integral to be evaluated is the volume of the triplets 
(x, y, 0) about the origin so that 

£ \ 2 / y £ \ 2 

- + (Jip m cos 0j + ( - + o-jp m sin 0± - j >p 2 . (12) 

As illustrated in figure El for different orientations of the vector connecting the two 
particles, this is a region bounded by four arc-circles, which we denote by 



y a ,r{x) = -2crp m sin - rl + r J 4p 2 - (x + 2ap m cos 0) 2 , (13) 
where a, r = ±1. 

As seen from figure O the area is connected for — 0t < < 0t, where 0t is the 
angle at which opposite arcs intersect, 

2 2 

T = arcsin Pm , - Pc . (14) 

lPm 

Beyond that value, the splits into two triangular These areas shrink to zero 

at the angle given by 



p m p c + 1/2^/ p 2 -pi 



ma- ' . - v ' ' Pm /, r \ 

= arccos r 2 . (15) 

P 
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Figure 5. Possible positions of the center of mass (x,y) for different values of </>, see 
equation (fTTj) , at a binary collision event. Here p = 13/25Z and p nl = 13/50L 



Let — 0m < < 4>m- We denote by x\ < %2 < £3 < £4 the four corners of the rectangular 
domain, 

x x = -x 4 = -2(p m cos0 - p c ), 
£2 = -x 3 = -2 a/p 2 - /4 sin 0. 

For cf) > 0T; the points at which the opposite arcs intersect are given by 

. ^ -/ 2 + 4/ Pm sin0 + 4(p 2 -p2^ 2 



Z 2 -4Zp m sin0 + 4p 2 n 

Combining equations ffl~3]) - ffTT|) together, we can make use of the symmetry 
and write the integral to be computed as 



(16) 



;i7) 



a(p, p r 



d(pdR 



A 1 (<f ) )d<f ) + / A 2 (0)d0 



where 
^i(0) 



£3 



j/ + i_i(a;)da; 



X"4 



y-i-i(x)dx 



X3 



Xl> 



y+i i+ i(a;)da; 



X4 



X2 



which is the area bounded by the four arcs y aT , and 



A 2 {<f>) 



[y-i,+i(a:) - y+i,_i(a:)]da:, 



y-i,+i(x)dx, (19) 



(20) 



is the area of the overlapping opposite arcs y~i,+i and which occurs when 

0t < < 0m [it gives a negative contribution to Al(0)]. 

The computation of these expressions is easily performed numerically, and the result 
shown in figure [6j 

Near the critical geometry, we expand the quantity (jl8|) in powers of the difference 

Pm - Pc, 



«(p>pm) = y^c»(p D 



(21) 



n=l 



The first two coefficients vanish, so that the leading term corresponds to n = 3. The 
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Figure 6. Area of binary collisions a(p, p m ) versus p m — p c for seven values of p 
ranging from 9/25 to 42/100 (I = 1/a/2). 



first few coefficients are derived in Appendix B and given by : 



C\ = 


c 2 = 0, 




128 Pc 


c 3 = 


3/ 2 ' 




256p c 2 


c 4 = 


3/ 4 ' 




16 / 1 


c 5 = 





(22) 

p c 16p 3 c 16p 4 



+ 5Z 2 + I 4 ° ^ 5/ 4 p c 
We further note that the computation of the two-cell 4- volume, |£ p , Pm (2)|, which, 
as noticed above, is approximated by the square of the single cell area \B P \ 2 , can 
be improved using equation fl2!l . Indeed it is easily seen that d|£ PiPm (2)|/dp m = 
— 2a(p, p m ), which implies 

oo 

Cn-1 



1^(2)1 = \B P ? ~ 2 —fa - P^- ( 23 ) 

71=4 

It is immediate to check that corrections |£ PiPm (2)| — \B P \ 2 = 0(p m — p c ) 4 . 
3. 3. Binary collision frequency 

Having computed the transition rates of the master equation with respect to the billiard 
geometry, we now turn to the computation of observables. The first quantity of interest, 
which can be readily computed from equation ([TO]) , is the binary collision frequency v^. 
This is an equilibrium quantity which, in the (global) microcanonical ensemble with 
energy E = ei + . . . + ejv, involves the two-particle energy distribution, 



,(«0,. . \ (N-l)(N-2) ^ e a + e b 



PP> (e a , e b ) = ^ ( 1 - ^ ) • (24) 



N-3 



and can be written as 



J de a de 6 d?7W / (e a ,e 6 |e a -?7,e 6 + ?7)P 2 (eq) (e a ,e; ) ). (25) 
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Taking the large N limit and letting E = Nk^T = N//3, we can write p!f Ci \e a) e b ) — 
(5 2 exp[— (5{e a + e b )][l + O^N^ 1 )]. Substituting this expression into the above equation 
and inserting the expression of W from equation (fTUj) . we obtain, after some calculations, 



k ^W'(I mR ) ll+ ° (N ~ 1)] - (26) 

This expression involves the geometric factor ( fl8l) in the first bracket. The leading term 
in the second bracket is the canonical expression of the binary collision frequency. The 
second term is a positive finite N correction, which is useful in that it shows that the 
binary collision frequency decreases to its asymptotic value as N — > oo. 

3.4- Rescaled master equation 

The binary collision frequency (|26p defines a natural dimensionless time scale for the 
stochastic process described by the master equation (|9]) with transition rates (fTDl) . The 
master equation can thus be converted to dimensionless form by rescaling the energies 
by a reference thermal energy and time by the corresponding asymptotic (N — ► oo) 
value of the binary collision frequency ( }26l) . 
Introducing the variables 

(27) 

h - £?• < 28 ' 

r =ubt, (29) 



equation (jUJ) becomes 

1 N 

^rp!i eq) ( e l ; • • • , e N, T) = - ^ 



2 

a,b=l 



dh 



X 



w(e a + h,e b - h\e a , e b ) p^ eq) (. . . , e a + h, . . . , e b - h, . . . , r) 
-w(e aj e b \e a - h,e b + h) p^ q \. . . ,e aj . . . ,e b ,. . . ,t) , (30) 
with the transition rates 

w{e a ,e b \e a - h, e b + h) = J / dx 1 \\dx 1 ±dx 2 \\dx 2 ±(x 1 \\ - x 2 \\) (31) 

II — "^2 1| 

x ^(e a — ^i|j — X\±) S(e b — x%h — x% ± ) S(h — x\» + x%») 

This master equation shows that all the properties of heat conduction are rescaled by 
the binary collision frequency and temperature in its limit of validity where the collision 
frequency vanishes. In particular, this shows that the coefficient of heat conduction 
is proportional to the binary collision frequency in this limit, as explained in the next 
subsection. 
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3.5. Thermal conductivity 

Starting from the master equation ([9]), we derive an equation for the evolution of the 
kinetic energy of each moving particle, (e a ) = k&Ta, which defines the local temperature, 
where (.) denotes an average with respect to the energy distributions. By the structure 
of equation ([9]), such an equation can be expressed in terms of the transfer of energy due 
to the binary collisions between neighbouring cells. The time evolution of the average 
local kinetic energy is given by 



This expresses the local conservation of energy. 

Over long time scales, the probability distribution becomes controlled by this local 
conservation of energy, the slowest variables being the local kinetic energies (e a ) or, 
equivalently, the local temperatures T a as defined above. This holds even though 
statistical correlations develop between the local energies in the probability distributions. 
These statistical correlations are well known for transport processes ruled by master 
equations such as Eq. ([9]) [19], [23] and are observed in the present system as well. 

To be specific, we consider a one-dimensional chain, extending along the x-axis, 
formed with a succession of pairs of rhombic billiard cells arranged in quincunx, similarly 
to the middle panel of figure HI except the vertical height is here only one unit of 
length. The unit of horizontal and vertical lengths is thus l\[2 and there are two cells 
per each unit of length. Similar results hold for different choices of geometry modulo 
straightforward adaptations. 

We imagine that the system is in a non-equilibrium state, with a small temperature 
difference 5T about an average temperature T between neighbouring cells, and consider 
the average heat transfered from cell a at inverse temperature (3 a = (3 + 5/3/2 to cell b 
at inverse temperature (3b = (3 — 5(3/2, 5(3 = — #T/(/c B T 2 ), both cells being assumed to 
be in thermal equilibrium at their respective temperatures. The statistical correlations 
we observe in the present system are of the order of 5(3 2 , as it is the case in other 
systems [19]. In the non-equilibrium state, these statistical correlations are controlled 
in the long-time limit by the local temperatures. Since the process is here ruled by a 
Markovian master equation in the limit of small binary collision frequency, we get the 
equation of heat for the temperature 




(32) 



b 



with the energy flux defined as 




(33) 



d t T(x,t) = d x [Kd x T(x,t)}, 

where the local temperature is here written 
i — 1, . . . , N/2, j — 0,1 [see equation ([I])]. 




(34) 
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According to the rescaling property of the master equation discussed in section 13. 4[ 
the heat conductivity is proportional to the binary collision frequency [Jj] : 



with a dimensionless constant A. 

An analytical estimation can be obtained by transforming the master equation into 
a hierarchy of equations for all the moments of the probability distribution: (e ), (e a e b ), 
(e a e&e c ),... The evolution equations of these moments are coupled. 

Truncating the hierarchy at the equations for the averages (e a ), we get the 
approximate heat conductivity : 



So that A = 1 in this approximation. The same result holds if we include the equations 
for the moments {e a e b ) with |a — b\ < 1. 

Though the approximate result fl36|) above does not rule out possible corrections to 
A = 1, we argue that A = 1 is an exact property of the master equation ([9]) in the infinite 
system limit. This claim is borne out by extensive studies of the stochastic process 
described by equation ([91) that will be reported in a separate publication [21]. The 
focus is here on the billiard systems whose conductivity may however bear corrections 
to this identity, due we believe to lack of sufficient separation of time scales between 
wall and binary collision events. In the following, the results of numerical computations 
are presented which support these claims. 

4. Numerical results 

The above formulae for the binary collision frequency, equation ( |26|) . and the thermal 
conductivity, equation (|36|) . together with the expressions of the geometric factors, 
equations (122!) and (|23j) . provide a detailed picture of the mechanism which governs the 
transport of heat in our model. Numerical computations of these quantities further add 
to this picture and provide strong evidence of the validity of our theoretical approach. 

4-1. Binary and wall collision frequencies 

For the sake of computing the binary collision frequency, we simulate the quasi-one 
dimensional channel of N cells with rhombic shapes and apply periodic boundary 
conditions at the horizontal ends of the channel. Thus each cell has two neighbouring 

I We notice that the heat capacity per particle is equal to cy = (\/N)dE/dT = k&, so that the thermal 
conductivity is also equal to the thermal diffusivity in units where fee = 1- 



(35) 




de a de fc dr? i](e b - e a )W(e a , e b \e a -n,e b + rj) exp[-/5(e a + e 6 )], 




(36) 
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cells, left and right, and interactions between any two neighbouring cells can occur 
through both top and bottom corners 0. 

In figure [7J we compare the computations of i>\, to u w for a system of N = 10 
cells at unit temperature. The parameters are taken to be p = 9/25 and p m = 
3/25, . . . , 17/50 by steps of 1/50. Both collision frequencies are computed in the units 
of the microcanonical average velocity, Un = 2 N \/ 'N/2(N — 1)\/(2N — 1)!!, which, as 
N — > oo, converges to the canonical average velocity Uqo = The wall collision 

frequency is compared to the collision frequency of the isolated cells u c , itself measured 
in the units of the single particle velocity. As expected, z/ b /z/ w <C 1 and z/ w ~ u c for 
p m > p c . The crossover z/b ~ v w occurs at p m ~ 11/50 for this value of p. 

2 
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Figure 7. Wall and binary collision frequencies v w and h>h versus p m — p c in a one- 
dimensional channel of N — 10 rhombic cells with p — 9/25 and lattice spacing 
I = 1/V2. 
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4-2. Thermal conductivity 

4-2.1. Heat flux. The thermal conductivity can be obtained by computing the heat 
flux in a nonequilibrium stationary state. Such stationary states occur when the two 
ends of the channel are put in contact with heat baths at separate temperatures, say 
T_ and T + , T_ < T + . 

Let a system of iV cells be in contact with two thermostated cells at respective 
temperatures T±, and let these cell indices be n = ±(N + l)/2 (we take N odd for the 
sake of definiteness) . Provided the difference between the bath temperatures is small, 
T + — T_ <^ (T + + T_)/2, a linear gradient of temperature establishes through the system, 
with local temperatures 

I TL 

T n = -(T+ + T_) + WTT (T + - T_). (37) 

% We mention that this set-up allows for re-collision between two particles (under stringent conditions), 
due to the vertical periodic boundary conditions. This conflicts with the assumptions of [15j . but does 
not seem to affect the results as far as numerics are concerned. 
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Under these conditions, the nonequilibrium stationary state is expected to be locally 
well approximated by a canonical equilibrium at temperature T n . 

This can be checked numerically. In fact the local thermal equilibrium is verified (by 
comparing the moments of e n to their Gaussian expectation values) under the weaker 
property of small local temperature gradients, i. e. (T+ — TJ)/N <C (T+ + T_)/2, for 
which the temperature profile is generally not linear since the thermal conductivity 
depends on the temperature. Indeed, since k oc T 1 / 2 , we expect in that case, according 
to Fourier's law, the profile 

2/3 

(38) 

The thermal conductivity can therefore be computed from the heat exchanges of 
the chain in contact with the two cells at the ends of the chain, respectively thermalized 
at temperatures T ± 5T/2, 5T <C T. The thermalization of the end cells is achieved by 
randomizing the velocities of the two particles at every collision they make with their 
cells walls, according to the usual thermalization procedure of particles colliding with 
thermalized walls [28] . 

First, we consider a chain containing a single cell in order to test the validity of 
the master equation. In this case, the procedure amounts to simulating a single particle 
confined to its cell and performing random collisions with stochastic particles which 
penetrate the cell corners according to the statistics of binary collisions. For a chain 
with a single cell, the heat conductance is given by Eq. ( 1331) with A = 1, as there are 
no correlations with the stochastic particles. 

Figure M shows the results of the computations of the heat conductance with this 
method and provides a comparison with the binary collision frequency on the one 
hand (left panel), as well as with the results of our kinetic theory predictions on the 
other hand (right panel). 

The agreement between the data and equations (1231) . (133|) . (1331) and the computation 
of the integrals (1X51) -( 1231) . especially as p m — > p c , demonstrates the validity of the 
stochastic description of the billiard system, equation (jUJ). 

Next, we increase the size of the chain in order to reach the thermal conductivity 
in the limit of an arbitrarily large chain. The results are that statistical correlations 
appear between the kinetic energies along the chain. As we show below, their influence 
on the computed value of the conductivity diminishes as p m — ► p c . 

Fix T_ = 0.5 and T + = 1.5 to be the baths temperatures, and let the size of the 
system increase from N = 1 to iV = 20 as p m is progressively decreased from p m = 11/50 
to p m = 13/100, with fixed p = 9/25. As one can see from figure [8], this range of values 
of p m crosses over from a regime where the separation of time scales is not effective 

+ Here and in the sequel, the thermal conductance or conductivity are further divided by l 2 VT, where 
/ = l/v2 is the rhombic cell size, so as to eliminate its length and temperature dependences, thus 
defining the reduced thermal conductivity k* = n/(l 2 ^/T). In these expressions and from here on, we 
further set fee = 1. 
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Figure 8. Reduced thermal conductance k, computed from the heat exchange in a 
chain with a single particle with thermalized neighbours, and binary collision frequency 
j'b, as functions of p m — p c . (Left) ratio between k and v\,\ (Right) comparison with 
the results of section |3l The only relevant parameter is p = 9/25. For each value of 
p m , several temperature differences ST were taken, all giving consistent values of k. 
The solid line shows the result of kinetic theory (K.T.). 

to one where it appears to be and where the stochastic model should therefore be a 
reasonable approximation to the process of energy transport in the billiard. 

For all values of p m , we measured the temperature profile and heat fluxes throughout 
the system and inferred the value of the reduced thermal conductivity by linearly 
extrapolating the ratios between the average heat flux and local temperature gradient 
divided by the square root of the local temperature as functions of 1/N to the vertical 
axis intercept, corresponding to N = oo. Given the parameter values, every realisation 
was carried out over a time corresponding to 1,000 interactions between the system and 
baths and repeated over 10 4 realisations. For N up to 20, this time provided satisfying 
stationary statistics, with temperature profiles verifying equation (138|) and statistically 
constant heat fluxes. 

The results of the linear regression used to compute the value of n/vh for selective 
values of p m are shown in figures [9] and [T0l 

Our results thus make it plausible that the ratio between the thermal conductivity 
and binary collision frequency approaches unity as the parameter p m decreases towards 
the critical value p c . As will be show in a separate publication [24J, this is indeed 
a property of the stochastic model described by equation ([9]) and has been verified 
numerically by direct simulation of the master equation within an accuracy of 4 digits. 
As of the billiard system, it is unfortunately difficult to improve the results beyond 
those presented here as the CPU times necessary to either increase N or decrease p m 
quickly become prohibitive. Nevertheless, the data displayed in figures [9] and [ID] offer 
convincing evidence that the thermal conductivity is well approximated by the binary 
collision frequency so long as the separation of time scales between wall and binary 
collision events is effective. 
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Figure 9. Ratio between thermal conductivity and binary collision frequency, 
K/vb, extrapolated from the computation of the average ratio between heat current 
and local temperature gradients, divided by the local binary collision frequency, 
V^Ei^i+i/^i+i - Ti)]/v h (i). The system sizes are N = 1,2,5,10,15,20. The 
four pannels correspond to different values of p m = 13/100,4/25, 19/100, 11/50, with 
p = 9/25 and thus p c ~ 0.068. The red dots are the data points with corresponding 
error bars and the black solid line shows the result of a linear regression performed 
with data associated to systems of lengths N > 2. 
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Figure 10. Ratio between thermal conductivity and binary collision frequency, n/vh, 
computated as in figure O here collected for a larger set of values of p m . The horizontal 
axis shows the difference p m — Pc- The error bars are of the same order as those in 
figure [51 with deviations from unity of the data points of the same order as those in 
the latter figure. 
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4-2.2. Helfand moment. The computation of the thermal conductivity can also be 
performed in the global equilibrium microcanonical ensemble using the method of 
Helfand moments [25| 1261 127] . 

The Helfand moment has expression H(t) = J2 a x a (t)e a (t), where x a (t) denotes the 
horizontal position of particle a at time t and e a (£) = \ v a(t)\ 2 /2 its kinetic energy (the 
masses are taken to be unity). The computation of the time evolution of this quantity 
proceeds by discrete steps, integrating the Helfand moment from one collision event to 
the next, whether between a particle and the walls of its cell, or between two particles. 
Let {r n } ne z denote the times at successive collision events. In the absence of binary 
collisions, the energies are locally conserved and the Helfand moment changes according 
to H(r n ) = F(r„_i) + Y. a l x a( T n) ~ £ a (T„_i)] e a {T n -i) . If, on the other hand, a binary 
collision occurs between particles k and I, the Helfand moment changes by an additional 
term [xj. (r„) — x\ (r n )] [e* (r n + 0) — (r n — 0)] . Computing the time average of the squared 
Helfand moment, we obtain an expression of the thermal conductivity according to 

k = lim — I— l im J-([ZT(t b ) - H(r )] 2 ) (39) 

L^oo L^K-qI Y n->oo ZT n \ I 

where L = N/2 is the horizontal length of the system. 

Figure [11] shows the results of a computation of the thermal conductivity through 
equation (1391) for different system sizes. Though the actual values of k vary wildly with 
N, it is clear that a finite asymptotic value is reached for N ~ 10 2 . In this case, the 
constant of proportionality in equation ( 1331) takes the value A = 0.98 ± 0.08, close to 1. 
Similar results were obtained for other parameter values, and other cell geometries as 
well. 
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Figure 11. Reduced thermal conductivity, computed from the mean squared Helfand 
moment, versus 1/N. The parameters are p — 9/25, p m = 9/50. The system sizes 
vary from N = 4 to N = 100. The dashed line shows the binary collision frequency, 
zvb — 0.1225. The solid line shows a linear fit of the data, with y-intercept 0.12 ± 0.01, 
in agreement with the prediction (|35|) . 
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5. Lyapunov spectrum 

A key aspect of our model, which justifies the assumption of local equilibrium, is that 
it is strongly chaotic. This property can be illustrated through the computation of the 
Lyapunov spectrum and Kolmogorov-Sinai entropy in equilibrium conditions. 

As mentioned earlier, in the absence of interaction between the cells, p m < p c , The 
Lyapunov spectrum of a system of N cells has N positive and N negative Lyapunov 
exponents, which, if divided by the average speed of the particle to which they are 
attached, are all equal in absolute value. This reference value we denote by A + . The 
2N remaining Lyapunov exponents vanish. 

As we increase p m and let the particles interact, we expect that, in the regime 
< p m — Pc "C 1, where binary collision events are rare, the Lyapunov exponents will 
essentially be determined by A + multiplied by a factor which is specified by the particle 
velocities. The exchange of velocities thus produces an ordering of the exponents which 
can be computed as shown below. We note that the other half of the spectrum, which 
remains zero in this approximation, will only pick up positive values as a result of the 
interactions. 

Assume for the sake of the argument that N is large. The probability that a given 
particle with velocity v has exponent A = v\ + less than a value Aj = v^X + can be 
approximated by the probability that the particle velocity be less than Vi, which, if we 
assume a canonical form of the equilibrium distribution, is 

Prob(A < A*) = Prob(w < fj), 

rmvf/2 

= (3 deexp(— /3e), 

Jo 

= 1 - exp(-/5m^ 2 /2). (40) 

But this probability is simply (N — i + 1/2) /N. Therefore the half of the positive 
Lyapunov exponent spectrum, which is associated to the isolated motion of particles 
within their cells, becomes, in the presence of rare collision events, 

nr\ n i 1/2 

A i = A + ,/— hi-—- , i = l,...,N, (41) 
V mfj |_ % — 1/2 

with ordering Ai > A2 > . . . > \n- In particular, the largest exponent Ai grows 
like Vln N. The Kolmogorov-Sinai entropy on the other hand is extensive : hxs = 
N\ + ^ii/{2mf3). 

Refined expressions can be computed by taking the microcanonical distribution 
associated to a finite N. In particular, the expressions of the Lyapunov exponents 
become 



A; 




1 



1/2V 1 1/2 



N 



(42) 



We mention in passing that similar arguments are relevant and can be used to 
approximate the Lyapunov spectrum (actually half of it) of other models, such as a 
mixture of light and heavy particles 





Figure 12. (Top) Lyapunov exponents Xi versus i computed with a rhombic channel 
of size N = 10 cells and parameter p = 9/25. The different curves correspond to 
different values of p m = O.ff (bottom curve) to 0.34 (top curve) by steps of 0.01. 
The lines of stars are obtained for p m — 0.04 < p c , yielding the exponents A + and 
A_ associated to isolated cells, with all the particles at the same speed. The squares 
correspond to the first half of the spectrum as predicted by equation (|42|) . (Bottom) 
Xi versus p m — p c . The first of the two figures displays the first half of the positive 
part of the spectrum of exponents, Ai, . . . , Xn and compares them to the asymptotic 
estimate equation (|42"T) (straight lines) . The second plot shows the second half of the 
positive part of the spectrum, Ajv+i, • ■ ■ , A2jv-ij displaying their power-law scaling to 
zero as p m — > p c . 
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Figure 13. Extensivity of the Kolmogorov-Sinai entropy. The vertical axis shows 
the difference between the Kolmogorov Sinai entropy, here computed from the sum of 
the positive Lyapunov exponents for system sizes N = 3 (magenta circles) , 5 (blue up 
triangles), 10 (green squares), 15 (cyan diamonds), 20 (red down triangles), divided 
by the corresponding microcanonical average velocity, and the Lyapunov exponent of 
the isolated billiard cell at unit velocity. The curves fall nicely upon each other and 
converge to zero as p m — > p c . 



Let us again insist that equations ( I4TI) and f j42l account for only N of the 2N — 1 
positive Lyapunov exponents. N — 1 are vanishing within this approximation. For all 
these Lyapunov exponents, we expect the corrections to vanish with the binary collision 
frequency as we go to the critical geometry. 

Figures [12] and [13] show the results of a numerical computation of the whole 
spectrum of Lyapunov exponents and corresponding Kolmogorov-Sinai entropy for the 
one-dimensional channel of rhombic cells. The agreement with equation (H2l as p m — > p c 
is very good, at the exception perhaps of the last few among the first N exponents, whose 
convergence to the asymptotic value (H2l) appears to be slower. Interestingly, the largest 
exponents have a minimum which occurs at about the value of p m for which the binary 
and wall collision frequencies have ratio unity, see figure [7] Indeed, for larger radii p m , 
the spectrum is similar to that of a channel of hard discs (without obstacles) [29]. Note 
that the same holds for the ratio between the thermal conductivity and binary collision 
frequency, as seen from figure [8] Thus one can interpret the occurrence of a minimum 
of the largest exponent as evidence of a crossover from a near close-packing solid-like 
phase (p m < p) to a gaseous-like phase trapped in a rigid structure (p m > p c ). 

6. Conclusions 

To summarize, lattice billiards form the simplest class of Hamiltonian models for which 
one can observe normal transport of energy, consistent with Fourier's law. Geometric 
confinement restricts the transport properties of this system to heat conductivity alone, 
thereby avoiding the complications of coupling mass and heat transports, which are 
common to other many particle billiards. 
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The strong chaotic properties of the isolated billiard cells warrant, in a parametric 
regime where interactions among moving particles are seldom, the property of local 
equilibrium. This is to say, assuming that wall collision frequencies are order-of- 
magnitude larger than binary collision frequencies, that one is entitled to making 
a Markovian approximation according to which phase-space distributions are spread 
over individual cells. We are thus allowed to ignore the details of the distribution at 
the level of individual cells and coarse-grain the phase-space distributions to a many- 
particle energy distribution function, thereby going from the pseudo-Liouville equation, 
governing the microscopic statistical evolution, to the master equation, which accounts 
for the energy exchanges at a mesoscopic cell-size scale and local thermalization. The 
energy exchange process further drives the relaxation of the whole system to global 
equilibrium. 

This separation of scales, from the cell scale dynamics, corresponding to the 
microscopic level, to the energy exchange among neighbouring cells at the mesoscopic 
level, and to the relaxation of the system to thermal equilibrium at the macroscopic 
level, is characterized by three different rates. The process of relaxation to local 
equilibrium has a rate given by the wall collision frequency, much larger than the 
rate of binary collisions, which characterizes the rate of energy exchanges which 
accompany the relaxation to local thermal equilibrium, itself much larger than the 
hydrodynamic relaxation rate, given by the binary collision rate divided by the square 
of the macroscopic length of the system. 

On this basis, having reduced the deterministic dynamics of the many particles 
motions to a stochastic process of energy exchanges between neighboring cells, we are 
able to derive Fourier's law and the macroscopic heat equation. 

The energy transport master equation can be solved with the result, to be presented 
elsewhere [24] . that the binary collision frequency and heat conductivity are equal. 
Under the assumption that the wall and binary collision time scales of the billiard are 
well separated, the transposition of this result to the billiard dynamics is that : 

(i) The heat conductivity of the mechanical model is proportional to the binary 
collision frequency, i. e. the rate of collisions among neighbouring particles, 



evolution of probability densities is rigorously described by the master equation, 
and remains close to unity over a large range of parameter values for which we 
conclude the time scale separation is effective and the master equation therefore 
gives a good approximation to the energy transport process of the billiard ; 

(ii) The heat conductivity and the binary collision frequency both vanish in the limit 
of insulating system, p m — > p c , with (p m — p c ) 3 , 



where the coefficient C3 depends on the specific geometry of the binary collisions. 



- — = A, v h < u w , 
with a constant A that is exactly 1 at the critical geometry, p 



in 



p c , where the 
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Though both results are exact strictly speaking only in the limit p m — > p c , the first 
one, according to our numerical computations, is robust and holds throughout the range 
of parameters for which the binary collision frequency is much less than the wall collision 
frequency, <C v w . The deviations from A = 1 which we observed at intermediary 
values of p m , where the separation of time scales is less effective, are interpreted as 
actual deviations of the energy transport process of the billiard from that described by 
the master equation, where correlations between the motions of neighboring particles 
must be accounted for. 

Under the conditions of local equilibrium, the Lyapunov spectrum has a simple 
structure, half of it being determined according to random velocity distributions within 
the microcanonical ensemble, while the other half remains close to zero. The analytic 
expression of the Lyapunov spectrum that we obtained is thus exact at the critical 
geometry. The Kolmogorov-Sinai entropy is equal to the sum of the positive Lyapunov 
exponents and thus determined by half of them. It is extensive in the number of particles 
in the system, whereas the largest Lyapunov exponent grows like the square root of the 
logarithm of that number. As we mentioned earlier, we believe our method is relevant 
to the computation of the Lyapunov spectrum of other models of interacting particles 
[30J . The computation of the full spectrum, particularly regarding the effect of binary 
collisions on the exponents, remains an open problem. 
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Appendix A. Computation of the Kernel 

We provide in this appendix the explicit form of the transition rate W, given by equation 

We first substitute the two velocity integrals by two angle integrals, eliminating the 
two delta functions which involve only the local energies. 

dv a dv b e ab -v ab 5 [e a - ) 5 ( e b 
y/2 



mv„ \ „ / mvi 



J S (v- y [(eat • v a f - (e ab ■ v b ) 2 }^j 

/ d9 a d9 b (^cos9 a - ^cos9 b )5(r] - e a cos 2 9 a + e b cos 2 9 b ), (A.l) 
Jd+ 

where 9 a / b denote the angles of the velocity vectors v a / b with respect to the direction 
of the relative position vector joining particles a and b, e ab = (cos 0, sin 0), and the 
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angle integration is performed over the domain D + such that y/e^cos9 a > y^cos^. 

With the above expression flA.ll) . the explicit <p dependence has disappeared so 
that we have effectively decoupled the velocity integration from the integration over the 
direction of the relative position between the two colliding particles. We can further 
transform this expression in terms of Jacobian elliptic functions as follows. 

Let Xi = cos9i, i — a,b, in equation flA.lj) . which becomes 

[ ^ Xa [ dXb 9(y/e^x a - y/ebX b )(y/€^x a - y/e^x b )5(r]-e a x 2 a +e b xl) , (A. 2) 
m 5/2 J -i sjl-xl J -i \/\-x\ 

where 8(.) is the Heaviside step function. We thus have to perform the x a and x b 
integrations along the line defined by the argument of the delta function, 

V = £ a x 2 a - e b xl, (A.3) 

and that satisfies the condition 

Ve^x a > yJT b X b . (A.4) 

To carry out this computation, we have to consider the following alternatives : 



e a < e b , < r] < e a 



The solution of equation ( 1A.3I) which is compatible with equation (1A.4I) is 

Plugging this solution into equation (1A.2I) and setting the bounds of the Xb-integral 
to ±a/ (e a — n)/e b , the expression flA.2j) reduces to (omitting the prefactors) 



Xa=[ v_± I! dY\ (AJi) 



dx b - 2 7 == = ^=k( £ -^L), (A.6) 



where K denotes the Jabobian elliptic function of the first kind, 

»7r/2 



/•7T/2 

K{m) = / (1 — msin 6)- 1/2 d6 (m < 1). (A.7) 
Jo 



Thus the kernel is, in this case, 



W(e a , e b \e a - v ,e b + r,)= l Pm J —K (^—1 ) / <\<k\R. (A.8) 

K I^P,Pm( 2 )l \ me b \ e b 



[u) e a < e b , e a - e b < n < 



This case is similar to case (i), with equation (IA.5I) replaced by 

Xb= _(zi^Y 2 (A . 9) 



e b 

and —1 < x a < +1. The expression of the kernel corresponding to this case is 
therefore 



W Mea - + ,) = -^.J-^K (_^_ ) / WR. (A.10) 
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(iii) | e a < e b , -e b < 77 < e a - e b < 

This case is similar to case (ii), with — v(e& + v)/ e a < x a < +a/ (e& + v)/ e a- In this 
case, the expression of the kernel is given by 

W(e a , e b \e a - rj, e b + rf) = _ 2pm f f*!^ / d 0dR (A.ll) 

The cases with e a > e b are obtained from the cases above with the roles of a and b 
interchanged and i] — > —rj. 

Appendix B. Collision area near the critical geometry 

The reason why C\ and c 2 in equation (122]) vanish is that the angle difference is 0(p m —p c ) 
and the area Ai(cj)) = 0[(p m — p c ) 2 ]- These quantities are easily computed. 
Let p m = (1 + e)p c , e C 1. We have 

<f>T= 2 -^e-^ + 0(e^ (B.l) 

0M = ^e + ^ + O(e% (B.2) 

which indicates that the bounds of the angle integrals appearing in equation (fl8l) are 

0(e), with cf>M and 0t differing only to 0(e 2 ). 

The leading contribution to the integral a(p, p m ) therefore stems only from the 

integration of Ai(cj)), which we can compute explicitly by expanding p m about p c and 

taking into consideration that is 0(e). The result is 

"<I>m r<t>T 

A(0)d0~ / A(0)d0, 







64p 



s'\ (B.3) 



3/ 2 

which yields the leading coefficient c 3 in equation 

We can compute the coefficients of the next few powers in the expansion (pT!) in a 
similar fashion. First we notice that A 2 (4>) is 0(e 3 ) so that its integral between 0t and 
0m is 0(e 5 ). Therefore only the integral of ^4i(0) contributes to c 4 in equation (122]) . 

The computation of the next terms in the expansion is more involved since it 
requires the integration of A2 (</>), whose expression is : 



^2(0) = 8p arcsin 



p m Z sin 0- (p 2 m - pi) 



P 2 



1/2 



(B.4) 



- 4[p m Zsin0 - (p 2 m - p c 2 )] 1/2 (4p 2 n + I 2 - 4p m Zsin0) 1 / 2 . 

Expanding this expression for p m £-close to p c and e 2 -close to 0t, we get, to leading 
order, 

" A2(m ^w E , (B5) 
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Combining this expression with the 5th order contribution to the integral of Ai(4>), we 
obtain C5, equation fl22|) . We point out that this coefficient is actually much larger than 
C3 and C4, a reason being that negative powers of p c appear in its expression. The same 
holds of the next few coefficients. 
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